# -*- coding: utf-8 -*-
import numpy as np
from scipy import signal,misc,ndimage
from skimage import filters,feature,img_as_float,io,color,util 
from PIL import Image,ImageFilter
import matplotlib.pyplot as plt
from PIL.ImageFilter import (FIND_EDGES,EDGE_ENHANCE,EDGE_ENHANCE_MORE)
from skimage.transform import pyramid_gaussian,pyramid_laplacian,pyramid_reduce, pyramid_expand, resize


def plot_image(image,title=''):
    plt.title(title,size=8),plt.imshow(image)
    plt.axis('off')
    
#当像素值出现负数的时候
def any_neighbor_zero(img, i, j):
    for k in range(-1,2):
      for l in range(-1,2):
         if img[i+k, j+k] == 0:
            return True
    return False
#判断像素值大小
def zero_crossing(img):
  img[img > 0] = 1
  img[img < 0] = 0
  out_img = np.zeros(img.shape)
  for i in range(1,img.shape[0]-1):
    for j in range(1,img.shape[1]-1):
      if img[i,j] > 0 and any_neighbor_zero(img, i, j):
        out_img[i,j] = 255
  return out_img
#导数与梯度
#ker_x = [[-1, 1]]
#ker_y = [[-1], [1]]
#im = color.rgb2gray(io.imread('F:\\pythonlx\\图像实战练习\\images\\chairs.png'))
#加载老虎图片
#im=color.rgb2gray(io.imread('F:\\pythonlx\\图像实战练习\\images\\tiger.jpg'))
#print(np.max(im))
#添加v随机噪声,#噪声对梯度的影响
#sigma = 1
#sign = np.random.random(im.shape)
#sign[sign <= 0.5] = -1
#sign[sign > 0.5] = 1

#im=im+util.random_noise(im, var=sigma**2)
##im = filters.gaussian(im, sigma=0.25)
#im_x = signal.convolve2d(im, ker_x, mode='same')
#im_y = signal.convolve2d(im, ker_y, mode='same')
#
#im_mag = np.sqrt(im_x**2 + im_y**2)
#im_ang = np.arctan(im_y/im_x)
#plt.gray()
#plt.figure(figsize=(10,10))
#plt.subplot(231)
#plt.imshow(im)
#plt.title('original', size=12)
#plt.axis('off')
#plt.subplot(232)
#plt.imshow(im_x)
#plt.title('grad_x', size=12)
#plt.axis('off')
#plt.subplot(233)
#plt.imshow(im_y)
#plt.title('grad_y', size=12)
#plt.axis('off')
#plt.subplot(234)
#plt.imshow(im_mag)
#plt.title('||grad||', size=12)
#plt.axis('off')
#plt.subplot(235)
#plt.imshow(im_ang)
#plt.title(r'$\theta$', size=12)
#plt.axis('off')
#plt.subplot(236)
##绘制边缘强度
#plt.plot(range(im.shape[1]), im[0,:], 'b-', label=r'$f(x,y)|_{x=0}$', linewidth=5)
#plt.plot(range(im.shape[1]), im_x[0,:], 'r-', label=r'$grad_x (f(x,y))|_{x=0}$')
#plt.title(r'$grad_x (f(x,y))|_{x=0}$', size=12)

#在同一图像中显示大小和梯度
#im=np.zeros((im.shape[0],im.shape[1],3))
#im[...,0]=im_mag*np.sin(im_ang)
#im[...,1]=im_mag*np.cos(im_ang)
#plt.imshow(im)
#plt.title(r'||grad||+$\theta$', size=30)
#plt.legend(prop={'size': 12})
#plt.show()

#拉普拉算子----用于边缘检测
#中心系数为正的拉普拉斯算子
#ker_laplacian = [[0,-1,0],[-1, 4, -1],[0,-1,0]]
##加载椅子图片
#im = color.rgb2gray(io.imread('F:\\pythonlx\\图像实战练习\\images\\chairs.png'))
#im1 = np.clip(signal.convolve2d(im, ker_laplacian, mode='same'),0,1)
#plt.gray()
#plt.figure(figsize=(11,10))
#plt.subplot(121)
#plt.imshow(im)
#plt.title('original', size=20)
#plt.axis('off')
#plt.subplot(122)
#plt.imshow(im1)
#plt.title('laplacian convolved', size=20)
#plt.axis('off')
#plt.show()


#锐化和反锐化掩模
#1拉普拉斯滤波器锐化图像
#输入图像
#im=color.rgb2gray(io.imread('F:\\pythonlx\\图像实战练习\\images\\man.png'))
##应用
#im1=np.clip(filters.laplace(im)+im,0,1)
#plt.figure(figsize=(10,10))
#plt.subplot(121),plot_image(im,'original image')
#plt.subplot(122),plot_image(im1,'sharpened image')
#plt.tight_layout()
#plt.show()

#反锐化掩模----锐化图像=原始图像+（原始图像-模糊蚃）*总数
#使用Scipy的ndimage模块实现反锐化掩模
#定义rgb每个通道在[0,1]
#def rgb2gray(im):
#    return np.clip(0.2989 * im[...,0] + 0.5870 * im[...,1] + 0.1140 * im[...,2], 0, 1)


#im = rgb2gray(img_as_float(io.imread('F:\\pythonlx\\图像实战练习\\images\\man.png')))
#
##模糊图像
#im_blurred = ndimage.gaussian_filter(im, 5)
#plt.imshow(im_blurred)
##细节图像
#im_detail = np.clip(im - im_blurred, 0, 1)
#plt.imshow(im_detail)
##画图
#fig, axes = plt.subplots(nrows=2, ncols=3, sharex=True, sharey=True, figsize=(10, 10))
#axes = axes.ravel()
#axes[0].imshow(im)
#axes[0].set_title('Original image', size=15)
#axes[1].imshow(im_blurred)
#axes[1].set_title('Blurred image, sigma=5', size=15)
#axes[2].imshow(im_detail)
#axes[2].set_title('Detail image', size=15)
#alpha = [1, 5, 10]
#for i in range(3):
#    #锐化图像
#    im_sharp = np.clip(im + alpha[i]*im_detail, 0, 1)
#    axes[3+i].imshow(im_sharp)
#    axes[3+i].set_title('Sharpened image, alpha=' + str(alpha[i]), size=15)
#
#for ax in axes:
#    ax.axis('off')
#fig.tight_layout()
#plt.show()
#
#for i in range(1,25):
#    im_sharp = np.clip(im + i*im_detail, 0, 1)
#    plt.figure(figsize=(10,10))
#    plt.imshow(im_sharp)
#    plt.title('Sharpened image, alpha=' + str(i), size=30)
#    plt.axis('off')
##    plt.savefig('F:\\pythonlx\\图像实战练习\\images\\man.png、、sharpened_me_' + str(i).zfill(3) + '.jpg', bbox_inches='tight')
#    plt.close()

#索贝尔边缘检测算法
    #转为灰度图片
#im = color.rgb2gray(io.imread('F:\\pythonlx\\图像实战练习\\images\\house2.png')) 
#plt.gray()
#plt.figure(figsize=(10,11))
#plt.subplot(2,2,1)
#plot_image(im,'origial')
#plt.subplot(2,2,2)
#edges_x = filters.sobel_h(im)
#plot_image(edges_x,'sobel_x') 
#plt.subplot(2,2,3)
#edges_y = filters.sobel_v(im)
#plot_image(edges_y,'sobel_y')
#plt.subplot(2,2,4)
#edges = filters.sobel(im)
#plot_image(edges,'sobel')
#plt.subplots_adjust(wspace=0.1, hspace=0.1)
#plt.show()

#应用不同边缘检测滤波器
#plt.gray()
#plt.figure(figsize=(10,11))
#plt.subplot(3,2,1),plot_image(im,'origial')
#edges=filters.roberts(im)#罗伯特
#plt.subplot(3,2,2),plot_image(edges,'robots')
#
#edges = filters.scharr(im)
#plt.subplot(3,2,3),plot_image(edges,'scharr') 
#
#edges= filters.sobel(im)
#plt.subplot(3,2,4),plot_image(edges,'sobel')
#
#edges = filters.prewitt(im)
#plt.subplot(3,2,5),plot_image(edges,'prewitt')
#
#edges =np.clip(filters.laplace(im),0,1)
#plt.subplot(3,2,5),plot_image(edges,'laplace')
#
#plt.subplots_adjust(wspace=0.1, hspace=0.1)
#plt.show()

#Canny边缘检测器
#im =  ndimage .gaussian_filter(im, 4)
#im += 0.05 * np.random.random(im.shape)
#
#edges1 = feature.canny(im)
#edges2 = feature.canny(im, sigma=3)
#
##绘制图像
#fig, (ax1, ax2, ax3) = plt.subplots(nrows=1, ncols=3, figsize=(10, 12),
#                                    sharex=True, sharey=True)
#
#ax1.imshow(im, cmap=plt.cm.gray)
#ax1.axis('off')
#ax1.set_title('noisy image', fontsize=15)
#
#ax2.imshow(edges1, cmap=plt.cm.gray)
#ax2.axis('off')
#ax2.set_title('Canny filter, $\sigma=1$', fontsize=15)
#
#ax3.imshow(edges2, cmap=plt.cm.gray)
#ax3.axis('off')
#ax3.set_title('Canny filter, $\sigma=3$', fontsize=15)
#
#fig.tight_layout()
#
#plt.show()

#LoG滤波器和DoG滤波器

#def plot_kernel(kernel,s,name):
#    plt.imshow(kernel,cmap='YlOrRd')#cmap='jet')#'gray_r')
#    ax=plt.gca()
#    ax.set_xticks(np.arange(-0.5,kernel.shape[0],2.5))
#    ax.set_yticks(np.arange(-0.5,kernel.shape[1],2.5))
#    plt.colorbar()
#
#def LOG(k=12,s=3):
#    n=2*k+1
#    kernel=np.zeros((n,n))
#    for i in range(n):
#        for j in range(n):
#            kernel[i,j]=-(1-((i-k)**2)/(2.*s**2))*np.exp(-((i-k)**2+(j-k)**2)/(2.*s**2))/(np.pi*s**4)
#    kernel=np.round(kernel/np.sqrt((kernel**2).sum()),3)
#    return kernel
#
#def DOG(k=12,s=3):
#    n=2*k+1
#    sl,s2=s*np.sqrt(2),s/np.sqrt(2)
#    kernel=np.zeros((n,n))
#    for i in range(n):
#        for j in range(n):
#           kernel[i,j]=np.exp(-((i-k)**2+(j-k)**2)/(2.*sl**2))/(2*np.pi*sl**2)-np.exp(-((i-k)**2+(j-k)**2)/(2*np.pi*s2**2))
#    kernel=np.round(kernel/np.sqrt((kernel**2).sum()),3)
#    return kernel
#
#s=3
#img=color.rgb2gray(misc.imread('F:\\pythonlx\\图像实战练习\\images\\flowers2.png'))
#kernel=LOG()
#outimg=signal.convolve2d(img,kernel)
#plt.figure(figsize=(10,10))
#plt.subplot(221),plt.title('LOG kernel',size=8)
#
#plot_kernel(kernel,s,'LOG')
#plt.subplot(222),plt.title('output image',size=8)
#plt.imshow(np.clip(outimg,0,1),cmap='gray')
#
#kernel=DOG()
#outimg=signal.convolve2d(img,kernel)
#plt.figure(figsize=(10,10))
#plt.subplot(223),plt.title('DOG kernel',size=8),plot_kernel(kernel,s,'DOG')
#plt.subplot(224),plt.title('output image',size=8)
#plt.imshow(np.clip(outimg,0,1),cmap='gray')
#plt.show()



#img = color.rgb2gray(misc.imread('F:\\pythonlx\\图像实战练习\\images\\zebra.png'))
#Scipy ndimage 模块的LoG滤波器
#fig=plt.figure(figsize=(10,9))
#plt.gray()
#for sigma in range(1,10):
#    plt.subplot(3,3,sigma)
#    img_log=ndimage.gaussian_laplace(img,sigma=sigma)
#    plt.imshow(np.clip(img_log,0,1)),plt.axis('off')
#    plt.title('LOG with sigma='+str(sigma),size=10)
#plt.show()


#采用零交叉点计算的Marr-Hildreth边缘检测算法
#fig = plt.figure(figsize=(10,10))
#plt.gray() 
#for sigma in range(2,10, 2):
#    plt.subplot(2,2,sigma/2)
#    result = ndimage.gaussian_laplace(img, sigma=sigma)
#    result = zero_crossing(result)
#    plt.imshow(result)
#    plt.axis('off')
#    plt.title('LoG with zero-crossing, sigma=' + str(sigma), size=10)
#
#plt.tight_layout()
#plt.show()


#基于PIL发现和增强边缘
#加载图片
#im=Image.open('F:\\pythonlx\\图像实战练习\\images\\panda.png')
#plt.figure(figsize=(10,10))
#plt.subplot(2,2,1)
#plot_image(im,'Original(UMBC library)')
#i=2
#for f in (FIND_EDGES,EDGE_ENHANCE,EDGE_ENHANCE_MORE):
#    plt.subplot(2,2,i)
#    im1=im.filter(f)
#    plot_image(im1,str(f))
#    i+=1
#plt.show()

#图像金字塔---融合图像
#1.scikit-image transform pygramid 模块的高斯金字塔


#
#image = io.imread('F:\\pythonlx\\图像实战练习\\images\\lena.jpg')
##image = imread('../new images/apple.png')
#rows, cols, dim = image.shape
#高斯
##pyramid = tuple(pyramid_gaussian(image, downscale=2)) #, multichannel=True))
#拉普拉斯
#pyramid = tuple(pyramid_laplacian(image, downscale=2))
#plt.figure(figsize=(11,5))
#i, n = 1, len(pyramid)
#for p in pyramid:
#    plt.subplot(3,3,i)
#    plt.imshow(color.rgb2gray(p), cmap='gray')
##    plt.subplot(1,n,i)
##    plt.imshow(p)
#    plt.title(str(p.shape[0]) + 'x' + str(p.shape[1]))
#    plt.axis('off')
#    i += 1
##plt.suptitle('Gaussian Pyramid', size=10)
#
#plt.suptitle('Laplacian Pyramid', size=10)
#
#plt.show()
#
##composite_image = np.zeros((rows, cols + cols // 2, 3), dtype=np.double)
#composite_image = np.zeros((rows, cols + cols // 2), dtype=np.double)
##composite_image[:rows, :cols, :] = pyramid[0]
#composite_image[:rows, :cols] = color.rgb2gray(pyramid[0])
#
#i_row = 0
#for p in pyramid[1:]:
#    n_rows, n_cols = p.shape[:2]
##    composite_image[i_row:i_row + n_rows, cols:cols + n_cols] = p
#    composite_image[i_row:i_row + n_rows, cols:cols + n_cols] =color.rgb2gray(p)
#    i_row += n_rows
#
#fig, ax = plt.subplots(figsize=(8,8))
##ax.imshow(composite_image)
#ax.imshow(composite_image, cmap='gray')
#plt.show()

#构造高斯金字塔
'''
计算步骤;
1.从原始图像开始
2.迭代计算金字塔每一层的图像
3.首先平滑图像，然后对图像向下采样
4.在图像变得足够小的时候停止

'''
#算法实现程序
def get_gaussian_pyramid(image):
    rows, cols, dim = image.shape
    gaussian_pyramid = [image]
    while rows > 1 and cols > 1:
        #print(rows, cols)
        image = pyramid_reduce(image, downscale=2)
        gaussian_pyramid.append(image)
        #向下采样
        rows //= 2
        cols //= 2
    return gaussian_pyramid
'''
拉普拉斯计算步骤;
1.从高斯金字塔和最小图像开始
2.迭代计算当前层的图像和获得的图像，（先进行上采样，然后平滑高斯金字塔上一层的图像）之间的差值
3.图像大小和原始图像之间大小相等时停止
'''
#实现程序
def get_laplacian_pyramid(gaussian_pyramid):
    laplacian_pyramid = [gaussian_pyramid[len(gaussian_pyramid)-1]]
    for i in range(len(gaussian_pyramid)-2, -1, -1):
        image = gaussian_pyramid[i] - resize(pyramid_expand(gaussian_pyramid[i+1]), gaussian_pyramid[i].shape)
        #print(i, image.shape)
        laplacian_pyramid.append(np.copy(image))
    laplacian_pyramid = laplacian_pyramid[::-1]
    return laplacian_pyramid


#导入图像
#image=io.imread('F:\\pythonlx\\图像实战练习\\images\\bunny.png')
##调用
#gaussian_pyramid = get_gaussian_pyramid(image)
#laplacian_pyramid = get_laplacian_pyramid(gaussian_pyramid)
#    
#w, h = 20, 12
#for i in range(3):
#    plt.figure(figsize=(w,h))
#    p = gaussian_pyramid[i]
#    plt.imshow(p)
#    plt.title(str(p.shape[0]) + 'x' + str(p.shape[1]), size=20)
#    plt.axis('off')
#    w, h = w / 2, h / 2
#    plt.show()
#    
#w, h = 10, 6#为上面处理后高斯金字塔图像大小
#for i in range(1,4):
#    plt.figure(figsize=(w,h))
#    p = laplacian_pyramid[i]
#    plt.imshow(color.rgb2gray(p), cmap='gray')
#    plt.title(str(p.shape[0]) + 'x' + str(p.shape[1]), size=20)
#    plt.axis('off')
#    w, h = w / 2, h / 2
#    plt.show()

#仅仅通过拉普拉斯金字塔重建图像
def reconstruct_image_from_laplacian_pyramid(pyramid):
    i = len(pyramid) - 2
    prev = pyramid[i+1]
    plt.figure(figsize=(20,20))
    j = 1
    while i >= 0:
        prev = resize(pyramid_expand(prev, upscale=2), pyramid[i].shape)
        im = np.clip(pyramid[i] + prev,0,1)
        plt.subplot(3,3,j)
        plt.imshow(im)
        plt.title('Level=' + str(j) + ', ' + str(im.shape[0]) + 'x' + str(im.shape[1]), size=20)
#        plt.savefig('../my images/out/pyramid' + str(len(pyramid) - i).zfill(2) + '.png', bbox_inches='tight')
        prev = im
        i -= 1
        j += 1
    plt.subplot(3,3,j)
    plt.imshow(im)
    plt.title('Original image' + ' ' + str(im.shape[0]) + 'x' + str(im.shape[1]), size=20)
    plt.show()
    return im
#
#image = io.imread('F:\\pythonlx\\图像实战练习\\images\\parrot.jpg')[...,:3] / 255
#pyramid = get_laplacian_pyramid(get_gaussian_pyramid(image))
#im = reconstruct_image_from_laplacian_pyramid(pyramid)


#基于金字塔的融合
'''
三幅图像要大小相同
'''
def main(name):
    im=Image.open(name)
    im=im.resize((200,200),Image.NEAREST).convert('RGB')
    a=img_as_float(np.array(im))
    return a
A =main('F:\\pythonlx\\图像实战练习\\images\\fruits.png')

B =main('F:\\pythonlx\\图像实战练习\\images\\flowers2.png')
M = main('F:\\pythonlx\\图像实战练习\\images\\pepper.jpg')

#rows, cols, dim = A.shape
pyramidA = get_laplacian_pyramid(get_gaussian_pyramid(A))
pyramidB = get_laplacian_pyramid(get_gaussian_pyramid(B))
pyramidM = get_gaussian_pyramid(M)

#重建
pyramidC = []
for i in range(len(pyramidM)):
    im = pyramidM[i]*pyramidA[i] + (1-pyramidM[i])*pyramidB[i]
    #print(np.max(im), np.min(im), np.mean(im))
    pyramidC.append(im)

im = reconstruct_image_from_laplacian_pyramid(pyramidC)

plt.imshow(im)
plt.axis('off')
plt.show()
